In [1]:
import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import re
import sys
import cPickle
import glob
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
import common_functions
import half_nanoplate_functions as hnf

C:\Users\Scherer Lab E\Anaconda2\lib\site-packages\skimage\viewer\utils\core.py:10: UserWarning: Recommended matplotlib backend is `Agg` for full skimage.viewer functionality.
  warn("Recommended matplotlib backend is `Agg` for full "
In [3]:
cd "K:\Pat's_Projects\ParticleTrajectoryData"

K:\Pat's_Projects\ParticleTrajectoryData
In [4]:
store = pd.HDFStore('half_nanoplate_dynamics_processed.h5', mode='r')

Estimate External Force F

Using Stokes Law and the average velocity (of the arc) of the particle over the nanoplate to calculate the force of the trap on the particles.

In [5]:
'''Calculate the velocities over nanoplate and Stokes drag forces for each 
velocity. Calculate the mean and deviation in the forces for each L'''
keys = store.index['key']
results = pd.DataFrame()
results['L'] = store.index['L']
velocity_list = []
force_list = []
mean_force = []
std_force = []
for i,key in enumerate(keys):
    df = store.get(key).copy()
    velocities = hnf.calc_velocities_consec_frames_in_theta(df, only_over_plate=True)
    forces = hnf.calc_force_stokes_drag(velocities)
    velocity_list.append(velocities)
    force_list.append(forces)
    mean_force.append(forces.mean())
    std_force.append(forces.std())
In [6]:
'''Add columns for last frame on plate and first frame on glass in new 
DataFrames'''
kinetic_data = []
keys = store.index['key']
for key in keys:
    df = store.get(key).copy()
    kinetic_data.append(hnf.add_barrier_crossing_index(df))
In [7]:
'''Fitting functions for dwell times'''
single_exp = lambda x, p1,p2: p1*np.exp(-p2*x)
double_exp = lambda x1, x2, p1,p2,p3: p1*(np.exp(-p2*x1) - np.exp(-p3*x2))

Dwell Times using a defined cutoff (230 degrees)

Set the cutoff for the attempts histogram for all experiments to be the same at 230 degrees

In [8]:
all_expts_dwell_time = []
for num, v in enumerate(kinetic_data):
    # Apply selection criteria
    temp_df = v[v['over_plate']==True]
    temp_df = temp_df.drop_duplicates(['frame', 'track id'])
    temp_df = temp_df.query('230 < theta < 280')
    
    # Find the dwell time for each particle leaving the plate in the DataFrame
    target_indices = temp_df[temp_df['last_frame_plate'] == True].index
    one_expt_dwell_time = []
    for target in target_indices:
        dwell_time = hnf.trajectory_dwell_time_to_index(temp_df, target, max_frames_absent=5)
        one_expt_dwell_time.append(dwell_time)
    all_expts_dwell_time.append(one_expt_dwell_time)
In [9]:
cum_dwell_time = plt.figure(figsize=(8,8))
Ls = [1, 2, 3, 4, 5, 10]
theta_230_cutoff_dwell_time = pd.DataFrame()
theta_230_cutoff_dwell_time['L'] = results['L']
cum_dwell_times = {}
fit_dwell_times = {}
for num, v in enumerate(all_expts_dwell_time):
    theta_230_cutoff_dwell_time.loc[num, 'mean_dwell_time'] = np.asarray(v).mean()
    
    # Histogram
    bins = hnf.hist_bin_optimization(v)
    hist, bins = np.histogram(v, bins=bins)
    mid_bins = (bins[:-1] + bins[1:])/2.0
    d_t = 1 - np.cumsum(hist)/float(sum(hist))
    if num == 0:
        params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, p0=[0.5,0.0095], maxfev=2000)
    else:
        params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, maxfev=2000)
    fit_x = np.linspace(min(mid_bins)*0.90, max(mid_bins)*1.10, 300)
    fit_y = [single_exp(xp, params[0], params[1]) for xp in fit_x]
    theta_230_cutoff_dwell_time.loc[num, 'decay'] = params[1]
    
    # Put data in lists for later
    cum_dwell_times['cdf_l'+str(Ls[num])] = d_t
    cum_dwell_times['dwell_l'+str(Ls[num])] = mid_bins/90.0
    fit_dwell_times['cdf_fit_l'+str(Ls[num])] = fit_y
    fit_dwell_times['dwell_fit_l'+str(Ls[num])] = fit_x/90.0
    
    
    # Plot
    plt.subplot(3,2,num+1)
    plt.plot(mid_bins/90.0, d_t, 'o')
    plt.plot(fit_x/90.0, fit_y, 'k')
    plt.title(keys[num])
    plt.title('Cumulative Dwell Time L='+str(Ls[num]))
    plt.xlabel('Dwell Time (sec)')
    plt.ylabel('$1-\mathregular{CDF}$')    
plt.tight_layout()
plt.show()
# save_dir = "C:\Users\Scherer Lab E\Dropbox\SPIE Presentation"
# cum_dwell_time.savefig(save_dir+"\cum_dwell_time_dist.pdf", dpi=800)
# cum_dwell_time.savefig(save_dir+"\cum_dwell_time_dist.png", dpi=800)
theta_230_cutoff_dwell_time
Out[9]:
L mean_dwell_time decay
0 1 69.241379 0.011996
1 2 18.303371 0.066009
2 3 8.320611 0.173177
3 4 5.144781 0.330313
4 5 4.065789 0.520610
5 10 2.277098 0.945243
In [14]:
# Collect all model fit parameters in a dictionary
model_fit_params = {}

Bell Model (230 Degree Cutoff)

In [28]:
k_f = theta_230_cutoff_dwell_time['decay']*90.0
mean_force_list = np.asarray(mean_force)
std_dev_force_list = np.asarray(std_force)
fit_params = np.polyfit(mean_force_list[:-1]*10**12, np.log(k_f)[:-1], 1)
print fit_params

print "k_o = "+str(np.exp(fit_params[1]))+" s^-1"
print "x_b = "+str(fit_params[0]*298*1.38*10**(-23)*10**12)+" meters"

k_o = np.exp(fit_params[1])
x_b = fit_params[0]*298*1.38*10**(-23)*10**12

model_fit_params['bell_k_o'] = k_o
model_fit_params['bell_x_b'] = x_b

bell_model_fig = plt.figure()
plt.plot(mean_force_list[:-1]*10**12, k_f[:-1], 'bo')
plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
fit_func = np.poly1d(fit_params)
plt.plot(np.linspace(0,0.25), np.exp(fit_func(np.linspace(0,0.25))), 'g')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('Dwell Time Decay Rate ($\mathregular{s}^{-1}$)')
plt.title('$k(F)$ vs Force')

chi2 = np.sum(((np.log(k_f)[:-1] - fit_func(mean_force_list[:-1]*10**12))**2)/(fit_params[0]*(std_dev_force_list[:-1]*10**12))**2)
print "Chi-squared fit: "+str(chi2)[:-1]
# save_dir = "C:\Users\Scherer Lab E\Dropbox\SPIE Presentation"
# bell_model_fig.savefig(save_dir+"\\bell_model_fit.pdf", dpi=800)
# bell_model_fig.savefig(save_dir+"\\bell_model_fit.png", dpi=800)

[ 29.2508446   -0.46045896]
k_o = 0.63099397843 s^-1
x_b = 1.20291173322e-07 meters
Chi-squared fit: 0.30078615259

Hummer Szabo Fit

There is a paper written by Hummer and Szabo that describes an extension of Bell's Model within the context of single-molecule pulling experiments. The reference is:

Intrinsic Rates and Activation Free Energies from Single-Molecule Pulling Experiments, DOI: 10.1103/PhysRev/Lett.96.108101

The main I am trying to fit to my data is eqn 3 in the text. The reason why is because in the paper in Fig 4 the force versus rate looks similar to what I see for my data. The equation being fit is (different from paper because it is missing some symbols):

$k(F) = k_o\left(1 - \frac{\nu F x^\ddagger}{\Delta G^\ddagger}\right)^{\frac{1}{\nu}-1} e^{\beta \Delta G^\ddagger[1 - (1 - \nu F x^\ddagger / \Delta G^\ddagger)^{1/\nu}]}$

Where $k(F)$ is the rate as a function of force ($F$), $k_o$ is the rate at $F=0$, $x^\ddagger$ is the distances from the lowest point in the well to the transition state ($\Delta G^\ddagger$). $\nu$ describes the nature of the barrier. $\nu=1/2$ corresponds to the cusp free-energy surface while $\nu = 2/3$ corresponds to a linear-cubic free energy surface. At $\nu=1$ the equation reduces to the Bell Model. I try fitting this equation using $k_o, \Delta G^\ddagger,$ and $x^\ddagger$ as fitting parameters for different values of $\nu$.

Fit when $\nu$ = 1/2

In [29]:
def hum_szabo_unified_theory(force, k0, xdd, gdd):
    first_multiple = k0 * (1.0 - (0.5*force*xdd/gdd))**(1.0/0.5 - 1.0)
    exponent = np.exp((gdd/kbt)*(1.0 - (1.0 - (0.5*force*xdd/gdd))**(1.0/0.5)))
    return first_multiple * exponent
In [30]:
kbt = 1.3806*10**(-23)*298
k_f = theta_230_cutoff_dwell_time['decay']*90.0
input_force = mean_force_list[:-1]
print k_f
fit_params, cov = scipy.optimize.curve_fit(hum_szabo_unified_theory, input_force, k_f[:-1], maxfev=2000000000, p0=[3*10**-4, 300*10**-9, 12*kbt])

print "k_o = "+str(fit_params[0])+" s^-1"
print "xdd = "+str(fit_params[1])+" meters"
print "del_G = "+str(fit_params[2]/kbt)+" kt"
print fit_params
print mean_force_list
plt.plot(mean_force_list[:-1]*10**12, k_f[:-1], 'bo')
plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
fit_x = np.linspace(0,0.25*10**-12)
fit_y = hum_szabo_unified_theory(fit_x, *fit_params)
plt.plot(fit_x*10**12, fit_y, 'g')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('Dwell Time Decay Rate')
plt.title('$k(F)$ vs Force')
plt.xlim(0,0.25)

theo_y = hum_szabo_unified_theory(input_force, *fit_params)
theo_sigma = hum_szabo_unified_theory(std_dev_force_list[:-1], *fit_params)
chi2 = np.sum(((k_f[:-1] - theo_y)**2)/(theo_sigma)**2)
print "Reduced Chi-squared fit: "+str(chi2)[:-1]

model_fit_params['dhs_n0.5_k_o'] = fit_params[0]
model_fit_params['dhs_n0.5_xdd'] = fit_params[1]
model_fit_params['dhs_n0.5_del_g'] = fit_params[2]/kbt

0     1.079607
1     5.940799
2    15.585965
3    29.728175
4    46.854939
5    85.071889
Name: decay, dtype: float64
k_o = 0.220516212774 s^-1
xdd = 2.51103469525e-07 meters
del_G = 7.26063657253 kt
[  2.20516213e-01   2.51103470e-07   2.98716239e-20]
[  3.14026011e-14   6.61943146e-14   1.01336498e-13   1.25411851e-13
   1.59273529e-13   1.81012608e-13]
Reduced Chi-squared fit: 0.37768193182

Fit when $\nu$ = 2/3

In [31]:
def hum_szabo_unified_theory(force, k0, xdd, gdd):
    first_multiple = k0 * (1.0 - ((2/3.0)*force*xdd/gdd))**(1.0/(2/3.0) - 1.0)
    exponent = np.exp((gdd/kbt)*(1.0 - (1.0 - ((2/3.0)*force*xdd/gdd))**(1.0/(2/3.0))))
    return first_multiple * exponent
In [32]:
k_f = theta_230_cutoff_dwell_time['decay']*90.0
input_force = mean_force_list[:-1]
print k_f
fit_params, cov = scipy.optimize.curve_fit(hum_szabo_unified_theory, input_force, k_f[:-1], maxfev=2000000000, p0=[1e-01,  1.28915220e-07, 3.05213988e-20])

print "k_o = "+str(fit_params[0])+" s^-1"
print "xdd = "+str(fit_params[1])+" meters"
print "del_G = "+str(fit_params[2]/kbt)+" kt"
print fit_params
hum_sza_fit = plt.figure()
plt.plot(mean_force_list[:-1]*10**12, k_f[:-1], 'bo')
plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
fit_x = np.linspace(0,0.25*10**-12)
#fit_params = [1e-01,  1.28915220e-07, 3.05213988e-20]
fit_y = hum_szabo_unified_theory(fit_x, *fit_params)
plt.plot(fit_x*10**12, fit_y, 'k')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('Dwell Time Decay Rate ($\mathregular{s}^{-1}$)')
plt.title('$k(F)$ vs Force')
plt.xlim(0, 0.25)

theo_y = hum_szabo_unified_theory(input_force, *fit_params)
theo_sigma = hum_szabo_unified_theory(std_dev_force_list[:-1], *fit_params)
chi2 = np.sum(((k_f[:-1] - theo_y)**2)/(theo_sigma)**2)
print "Reduced Chi-squared fit: "+str(chi2)[:-1]
# save_dir = "C:\Users\Scherer Lab E\Dropbox\SPIE Presentation"
# hum_sza_fit.savefig(save_dir+"\hum_sza_fit.pdf", dpi=800)
# hum_sza_fit.savefig(save_dir+"\hum_sza_fit.png", dpi=800)

model_fit_params['dhs_n0.66_k_o'] = fit_params[0]
model_fit_params['dhs_n0.66_xdd'] = fit_params[1]
model_fit_params['dhs_n0.66_del_g'] = fit_params[2]/kbt

0     1.079607
1     5.940799
2    15.585965
3    29.728175
4    46.854939
5    85.071889
Name: decay, dtype: float64
k_o = 0.331106530809 s^-1
xdd = 2.03640794117e-07 meters
del_G = 6.27201003815 kt
[  3.31106531e-01   2.03640794e-07   2.58042284e-20]
Reduced Chi-squared fit: 0.23356828193

The reason for the cutoff of the curve at .37 pN is because the term in the exponent $(1 - \nu F x^\ddagger / \Delta G^\ddagger)^{1/\nu}$ becomes imaginary. For certain values of $F$ the term $1 - \nu F x^\ddagger / \Delta G^\ddagger$ becomes negative but it is being raised to an exponent of 3/2 (when $\nu$=2/3). When $\nu$=1/2 this problem doesn't occur because the term will always be squared.

Fit when $\nu$ = 1

I couldn't get the fit for the Bell Model to work with $\nu$ = 1

In [33]:
def hum_szabo_unified_theory(force, k0, xdd, gdd):
    first_multiple = k0 * (1.0 - (1.0*force*xdd/gdd))**(1.0 - 1.0/1.0)
    exponent = np.exp((gdd/kbt)*(1.0 - (1.0 - (1.0*force*xdd/gdd))**(1.0/1.0)))
    return first_multiple * exponent
In [34]:
k_f = theta_230_cutoff_dwell_time['decay']*90.0
input_force = mean_force_list[:-1]
print k_f
fit_params, cov = scipy.optimize.curve_fit(hum_szabo_unified_theory, input_force, k_f[:-1], maxfev=2000000000, p0=[3*10**-4, 300*10**-9, 12*kbt])

print "k_o = "+str(fit_params[0])+" s^-1"
print "x_b = "+str(fit_params[1])+" meters"
print "del_G = "+str(fit_params[2])+" joules"
print fit_params
plt.plot(mean_force_list[:-1], k_f[:-1], 'bo')
plt.errorbar(mean_force_list[:-1], k_f[:-1], xerr=std_dev_force_list[:-1], linestyle='none')
fit_x = np.linspace(0,0.5*10**-12)
fit_y = hum_szabo_unified_theory(fit_x, *fit_params)
plt.plot(fit_x, fit_y, 'g')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('log(Dwell Time Decay Rate)')
plt.title('$log(k(F))$ vs Force')

0     1.079607
1     5.940799
2    15.585965
3    29.728175
4    46.854939
5    85.071889
Name: decay, dtype: float64
k_o = 0.00323820631674 s^-1
x_b = 2.48895082055e-07 meters
del_G = 3.17079036317e-12 joules
[  3.23820632e-03   2.48895082e-07   3.17079036e-12]
Out[34]:
<matplotlib.text.Text at 0x1f044198>

Compare Hummer Szabo $\nu$=1/2 to Bell Model

In [35]:
def hum_szabo_unified_theory(force, k0, xdd, gdd):
    first_multiple = k0 * (1.0 - ((1/2.0)*force*xdd/gdd))**(1.0/(1/2.0) - 1.0)
    exponent = np.exp((gdd/kbt)*(1.0 - (1.0 - ((1/2.0)*force*xdd/gdd))**(1.0/(1/2.0))))
    return first_multiple * exponent
In [36]:
bell_hum_szabo_graph = {}

k_f = theta_230_cutoff_dwell_time['decay']*90.0
mean_force_list = np.asarray(mean_force)
std_dev_force_list = np.asarray(std_force)
fit_params = np.polyfit(mean_force_list[:-1]*10**12, np.log(k_f)[:-1], 1)
print fit_params

print "k_o = "+str(np.exp(fit_params[1])*90.0)+" s^-1"
print "x_b = "+str(fit_params[0]*298*1.38*10**(-23)*10**12)+" meters"
bell_hum_fits = plt.figure()
plt.plot(mean_force_list[:-1]*10**12, k_f[:-1], 'bo')
plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
fit_func = np.poly1d(fit_params)
plt.plot(np.linspace(0,0.25), np.exp(fit_func(np.linspace(0,0.25))), 'g')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('Dwell Time Decay Rate')
plt.title('$k(F)$ vs Force')

bell_hum_szabo_graph['avg_force'] = mean_force_list[:-1]*10**12
bell_hum_szabo_graph['std_force'] = std_dev_force_list[:-1]*10**12
bell_hum_szabo_graph['k_f'] = k_f[:-1]
bell_hum_szabo_graph['bell_x'] = np.linspace(0,0.25)
bell_hum_szabo_graph['bell_y'] = np.exp(fit_func(np.linspace(0,0.25)))

chi2 = np.sum(((np.log(k_f)[:-1] - fit_func(mean_force_list[:-1]*10**12))**2)/(fit_params[0]*(std_dev_force_list[:-1]*10**12))**2)
print "Chi-squared fit: "+str(chi2)[:-1]

k_f = theta_230_cutoff_dwell_time['decay']*90.0
input_force = mean_force_list[:-1]
print k_f
fit_params, cov = scipy.optimize.curve_fit(hum_szabo_unified_theory, input_force, k_f[:-1], maxfev=2000000000, p0=[1e-01,  1.28915220e-07, 3.05213988e-20])

print "k_o = "+str(fit_params[0])+" s^-1"
print "xdd = "+str(fit_params[1])+" meters"
print "del_G = "+str(fit_params[2]/kbt)+" kt"
print fit_params
#plt.plot(mean_force_list[:-1]*10**12, k_f[:-1], 'bo')
#plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
fit_x = np.linspace(0,0.25*10**-12)
#fit_params = [1e-01,  1.28915220e-07, 3.05213988e-20]
fit_y = hum_szabo_unified_theory(fit_x, *fit_params)
plt.plot(fit_x*10**12, fit_y, 'k')
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('Dwell Time Decay Rate ($\mathregular{s}^{-1}$)')
plt.title('$k(F)$ vs Force')
plt.ylim(10**-1,10**2)

bell_hum_szabo_graph['hum_szabo_x'] = fit_x*10**12
bell_hum_szabo_graph['hum_szabo_y'] = fit_y

theo_y = hum_szabo_unified_theory(input_force, *fit_params)
theo_sigma = hum_szabo_unified_theory(std_dev_force_list[:-1], *fit_params)
chi2 = np.sum(((k_f[:-1] - theo_y)**2)/(theo_sigma)**2)
print "Reduced Chi-squared fit: "+str(chi2)[:-1]
# save_dir = "C:\Users\Scherer Lab E\Dropbox\SPIE Presentation"
# bell_hum_fits.savefig(save_dir+"\\bell_hum_fits.pdf", dpi=800)
# bell_hum_fits.savefig(save_dir+"\\bell_hum_fits.png", dpi=800)

[ 29.2508446   -0.46045896]
k_o = 56.7894580587 s^-1
x_b = 1.20291173322e-07 meters
Chi-squared fit: 0.30078615259
0     1.079607
1     5.940799
2    15.585965
3    29.728175
4    46.854939
5    85.071889
Name: decay, dtype: float64
k_o = 0.220516626627 s^-1
xdd = 2.51103355037e-07 meters
del_G = 7.26063474366 kt
[  2.20516627e-01   2.51103355e-07   2.98716163e-20]
Reduced Chi-squared fit: 0.37768084442
In [18]:
'''Save the k(F) vs Force for Bell and Hummer Szabo data for the paper'''
import cPickle

pickle_file = open("C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Data\\bell_hum_szabo_fit.pkl", 'w')
cPickle.dump(bell_hum_szabo_graph, pickle_file)
pickle_file.close()
In [37]:
'''Save the fit parameters from the Bell and Hummer Szabo models'''
import cPickle

pickle_file = open("C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Data\\bell_hum_szabo_model_fit_params.pkl", 'w')
cPickle.dump(model_fit_params, pickle_file)
pickle_file.close()
In [27]:
store.close()